この記事を読む前に：全体像を把握しよう
この記事では「地下水の流れ」と「水の化学反応」を同時にシミュレーションするプログラムの概念を説明します。
まず、登場する2つのツールの役割を確認しましょう。
| ツール | 役割 | たとえるなら |
|---|---|---|
| MODFLOW 6 | 地下水がどこをどの速さで流れるかを計算 | 川の流れをシミュレートする地図アプリ |
| PHREEQC | 水の中でどんな化学反応が起きるかを計算 | 水質を予測する化学実験室 |
この2つを Python でつなぐことを 「カップリング」 と呼びます。MODFLOW が「水の動き」を決め、その情報を PHREEQC に渡して「化学反応の結果」を受け取る、というキャッチボールを繰り返すイメージです。
地球化学とは、岩石・水・ガスの間で起きる化学反応を研究する学問です。たとえば「石灰岩が熱水に溶ける」「地下水に鉄が溶け出す」「温泉の成分が時間とともに変わる」といった現象を、数式とコンピュータで定量的に予測します。
前回のおさらい
第1回（モデル設計編）では、カラムモデルの構造と境界条件を設計しました。
カラムモデルとは何でしょうか？地下の岩石層を縦に切り取った、細長い試験管のようなモデルです。左端から水を注入して右端から出てくる水の化学組成がどう変わるかを追います。
注入口 排出口
↓ ↓
[水を押し込む] → [セル1][セル2]...[セル15] → [水が出てくる]
↑ここで化学反応が起きる↑今回はこのカラムモデルを実際に「動かす」コードを組み立てます。
全体の処理の流れ
プログラム全体は次の6ステップで動きます。それぞれ何をしているのかを先に把握しておくと、後の説明が理解しやすくなります。
Step 1: MODFLOW 6 で地下水の流れを計算する
↓
Step 2: 計算結果から「水がどのくらいの速さで流れているか」を取り出す
↓
Step 3: 「1回のループで何秒進むか」を自動で計算する
↓
Step 4: PHREEQC（化学反応エンジン）を準備する
↓
Step 5: カラム内の初期状態（水の組成・鉱物量）を設定する
↓
Step 6: ループを回して流れと反応を繰り返し計算するStep 1：必要なライブラリの準備
プログラムを書く前に、必要な「道具箱」を読み込みます。Python では、こうした道具箱を ライブラリ と呼びます。
import os
import numpy as np # 数値計算
import flopy # MODFLOW を Python から操作する
import pandas as pd # 表形式のデータを扱う
import matplotlib.pyplot as plt # グラフを描く
from phreeqpy.iphreeqc.phreeqc_dll import IPhreeqc # PHREEQCエンジンPHREEQC は本来、専用のソフトウェアとして動かすものです。phreeqpy はそれを Python から直接呼び出せるようにする「橋渡し役」のライブラリです。インストールは pip install phreeqpy の1行で完了します。
Step 2：モデルの数値設定
シミュレーションで使う数値をまとめて定義します。「なぜこの値なのか」も合わせて説明します。
カラムの形
| 変数 | 値 | 意味 |
|---|---|---|
L |
1.0 m | カラムの長さ |
N |
15 | カラムを15個のブロック（セル）に分割 |
dx |
1/15 ≈ 0.067 m | 1セルの幅 |
水の流れ方
| 変数 | 値 | 意味 |
|---|---|---|
K |
1.0×10⁻³ m/s | 透水係数（岩石の水の通しやすさ） |
phi |
0.3 | 間隙率（岩石の30%が空洞） |
Qin |
1.0×10⁻⁵ m³/s | 注入流量 |
- 透水係数：岩石が水をどれだけ通しやすいかを表す数値。砂や砂利は大きく、粘土は非常に小さい値になります。
- 間隙率：岩石の体積のうち、水が入れる空洞の割合です。
phi = 0.3は体積の30%が空洞であることを意味します。
Step 3：MODFLOW 6 で流れを計算する
パッケージの登録順序が重要
MODFLOW 6 でモデルを組み立てるとき、部品を追加する順番に決まりがあります。この順序を守らないとエラーになります。
正しい順序：
MFSimulation → TDIS → IMS → GWF → DIS → NPF → IC → WEL → CHD → OC
TDIS（時間設定）は必ず最初！MODFLOW 6 では「時間設定（TDIS）」が決まっていないと、「何ステップ分のデータを作ればいいか」がわからないため、後続のパッケージが正しく初期化できません。TDIS を先に作るのはそのためです。
各パッケージの役割
| パッケージ | 役割 |
|---|---|
| TDIS | シミュレーション期間と時間刻みを定義 |
| IMS | 方程式を解く「ソルバー」（計算エンジン） |
| GWF | 地下水流れモデル本体 |
| DIS | グリッド（格子）の形状を定義 |
| NPF | 透水係数など水の流れやすさを設定 |
| IC | シミュレーション開始時の水頭（水位）を設定 |
| WEL | 注入井（水を押し込む境界） |
| CHD | 一定水頭境界（出口の水圧を固定） |
| OC | 計算結果をどのファイルに保存するかを指定 |
組み立てたら sim.write_simulation() でファイルに書き出し、sim.run_simulation() で実際に計算を実行します。
run_simulation() を実行するには mf6（Windows では mf6.exe）という実行ファイルが必要です。MODFLOW の公式ページからダウンロードして、作業フォルダか PATH の通った場所に置いてください。
Step 4：計算結果から「流速」を取り出す
MODFLOW の計算が終わると、Cell Budget File（.cbb） という結果ファイルが生成されます。このファイルから「水がセルをどの速さで通過したか（Darcy 流速）」を取り出します。
# .cbb ファイルを開いて Darcy 流速を読み取る
cbb = flopy.utils.CellBudgetFile(cbb_path, precision="double")
spdis = cbb.get_data(text="DATA-SPDIS")[0]
v_darcy = float(np.mean(np.abs(spdis["qx"].flatten())))岩石全体の断面積を基準にした「見かけの流速」です。実際に水が移動する速度（間隙流速）は、これを間隙率で割って求めます。
間隙流速 = Darcy流速 ÷ 間隙率たとえば Darcy 流速が 3×10⁻⁴ m/s で間隙率が 0.3 なら、実際の水の速さは 1×10⁻³ m/s です。
Step 5：時間刻みを自動計算する（Courant 条件）
化学反応の計算を正確に行うには、「1ステップで水がちょうど1セル分だけ進む」ように時間刻みを設定する必要があります。これを Courant 数 = 1 の条件 と呼びます。
\[\Delta t = \frac{\Delta x}{v_{\text{pore}}}\]
- \(\Delta x\)：1セルの幅（例：0.067 m）
- \(v_{\text{pore}}\)：間隙流速（m/s）
v_pore = v_darcy / phi # 間隙流速
dt = dx / v_pore # Courant 数 = 1 を満たす時間刻み
n_steps = int(total_time / dt) # 総ステップ数水が1ステップで2セル以上飛び越えると、化学反応の計算に大きな誤差が生じます。「1ステップ＝1セル」にすることで、水の移動と化学反応をきれいに対応させることができます。この考え方を 数値安定性の条件 と言います。
Step 6：PHREEQC の準備
データベースの選択
PHREEQC は「熱力学データベース」と組み合わせて使います。どのデータベースを使うかで、計算できる温度・圧力の範囲が変わります。
| データベース | 対応温度 | 特徴 |
|---|---|---|
phreeqc.dat |
～25℃ | 標準的な地下水向け |
llnl.dat |
～300℃ | 深部熱水・地熱向け（今回使用） |
このシミュレーションは 150℃ の熱水を扱います。高温条件の熱力学データを持つ llnl.dat（Lawrence Livermore 国立研究所が開発）を選ぶ必要があります。
初期状態の設定
シミュレーション開始前に、カラム内の各セルに「初期の水の組成」と「鉱物の量」を設定します。
初期水質（全15セル共通）
| 項目 | 値 |
|---|---|
| 温度 | 150℃ |
| 圧力 | 100 atm |
| pH | 7.0 |
| Na・Cl | 各 1×10⁻³ mol/kgw |
初期鉱物（方解石 = 石灰石の主成分）
各セルに方解石を 0.5 mol 配置し、熱水に溶けていく様子をシミュレーションします。方解石の溶解速度は「速度論的モデル（KINETICS）」で計算します。
化学反応には「即座に平衡に達するもの」と「ゆっくり時間をかけて反応するもの」があります。速度論的モデル（Kinetics）は後者を扱い、「1秒間にどれだけ溶けるか」という反応速度を考慮します。実際の地下環境では多くの鉱物溶解は時間がかかるため、速度論モデルが欠かせません。
Step 7：カップリングループ（メインの繰り返し計算）
いよいよシミュレーションの核心です。以下の処理を n_steps 回繰り返します。
【1ステップの処理】
① SOLUTION 0（入口から入ってくる水の組成）を定義
② TRANSPORT（流れの条件）を設定
③ PHREEQC を実行 → 各セルの化学反応を計算
④ 結果（pH・Ca濃度・飽和指数）を記録
⑤ ①に戻るTRANSPORT ブロックの主な設定
| 設定項目 | 意味 |
|---|---|
-cells 15 |
セル数（カラムを15分割） |
-lengths 0.067 |
1セルの幅 [m] |
-time_step dt |
1ステップの時間 [s] |
-shifts 1 |
1ステップで水を1セル分移動させる |
-punch_cells 1-15 |
全セルの結果を出力する |
実装で特に間違えやすいポイントをまとめます。
① SOLUTION 0 を毎ステップ定義する SOLUTION 0 は「入口から入ってくる水」を表します。毎ステップ定義しないと、PHREEQC が上流境界を認識できずエラーになります。
② -lengths（複数形）を使う -length（単数形）は phreeqpy では無効です。必ず複数形の -lengths を使ってください。
③ -punch_cells 1-15（ハイフン区切り）にする -punch_cells 1 15 と書くと「セル1とセル15だけ」の意味になります。全セルを出力するには 1-15 のようにハイフンで範囲指定が必要です。
結果の見方
シミュレーションが終わると、以下のデータが results 辞書に蓄積されます。
| キー | 内容 |
|---|---|
time |
経過時間 [時間] |
pH_in / pH_out |
入口・出口の pH の時間変化 |
Ca_in / Ca_out |
Ca²⁺ 濃度の時間変化 [mol/kgw] |
SI_in / SI_out |
方解石の飽和指数（SI）の時間変化 |
spatial |
各時刻でのカラム内の空間プロファイル |
飽和指数（Saturation Index）は、ある鉱物が「溶けやすい状態か、沈殿しやすい状態か」を示す指標です。
- SI < 0：不飽和 → 鉱物が水に溶けていく
- SI = 0：平衡状態 → 溶解も沈殿もない
- SI > 0：過飽和 → 鉱物が沈殿しようとしている
今回のシミュレーションでは、SI が時間とともにどう変化するかを追跡します。
処理フローの全体像（まとめ）
main.py
│
├── 0. パラメータ定義（カラムサイズ・流量・温度など）
│
├── 1. MODFLOW 6 でモデルを構築・実行
│ └── 地下水の流れ場を計算する
│
├── 2. 計算結果から Darcy 流速を取得
│
├── 3. Courant 条件から時間刻み Δt を決める
│ └── Δt = セル幅 ÷ 間隙流速
│
├── 4. PHREEQC を初期化
│ ├── データベース読み込み（llnl.dat）
│ ├── 初期水質を設定（SOLUTION 1-15）
│ └── 鉱物量を設定（KINETICS 1-15）
│
└── 5. カップリングループ（n_steps 回繰り返す）
├── 入口水質を定義（SOLUTION 0）
├── 輸送・反応を計算（TRANSPORT + KINETICS）
├── 結果を取得（pH・Ca・SI）
└── 時系列データ・空間プロファイルを記録よくあるエラーと対処法
| エラーメッセージ | 原因 | 対処 |
|---|---|---|
FlopyException: Could not find nper |
WEL を TDIS より先に作成した | TDIS を最初に登録する |
Solution 0 is needed for transport |
SOLUTION 0 が未定義 | TRANSPORT ブロックに SOLUTION 0 を含める |
| 全セルの結果が出ない | -punch_cells 1 15 と書いてしまった |
-punch_cells 1-15 に変更（ハイフンで範囲指定） |
-length が無効 |
単数形を使った | -lengths（複数形）に変更 |
次回予告：解析・可視化編
第2回ではカップリングループを動かす仕組みを実装しました。次回（第3回）ではシミュレーション結果を使って以下を可視化します。
- 空間プロファイル：ある時刻でカラム全体の化学組成がどう分布しているか
- 飽和指数マップ：どこで方解石が溶けていて、どこで沈殿しているか
このシリーズの他の記事：
- #1 モデル設計編
- #2 実装編（カップリングの核心）（本記事）
- #3 解析・可視化編